Sierra Nevada Time Series

Data preparation script for Sierra Nevada time series analysis

Merge templates with 14C, C, and N data

Radiocarbon analyses for the 2001 samples were not run originally, but were completed on archived samples in 2020.

2001 mean radiocarbon profiles

Fig. 10. Mean profile \(\Delta\)¹⁴C for 2001 samples

Caption: Mean \(\Delta\)¹⁴C by depth for each site in 2001. Error bars show ±1 standard deviation, solid vertical line shows \(\Delta\)¹⁴C of the atmosphere in the year of sampling.

2009 radiocarbon profiles

Fig. 10. Profile \(\Delta\)¹⁴C for 2009 samples

Caption: Profile \(\Delta\)¹⁴C by depth for each site in 2009. Solid vertical line shows \(\Delta\)¹⁴C of the atmosphere in the year of sampling. Error bars not shown as only a single replicate profile was analyzed per site.

2019 mean radiocarbon profiles

Fig. 10. Mean profile \(\Delta\)¹⁴C for 2019 samples

Caption: Mean \(\Delta\)¹⁴C by depth for each site in 2019. Error bars show ±1 standard deviation, solid vertical line shows \(\Delta\)¹⁴C of the atmosphere in the year of sampling.

Change in \(\Delta\)¹⁴C over time between 2001 and 2019

Fig. 10. Mean profile \(\Delta\)¹⁴C for 2001 and 2019 samples

Caption: Mean \(\Delta\)¹⁴C by depth for each site in 2001 and 2019. Error bars show ±1 standard deviation. Vertical lines show \(\Delta\)¹⁴C of the atmosphere in 2001 (solid) and 2019 (dashed).

Incubation \(\Delta\)¹⁴C-CO₂ for 2019 bulk soils

Fig. 10. Rep 1 \(\Delta\)¹⁴C-CO₂ of 2019 bulk soil incubations

Caption: \(\Delta\)¹⁴CO₂ by depth for each site in 2019. (Still waiting for data from laboratory duplicates). Note that the 10-20 cm increment sample from the cold granite site appears to have been contaminated: \(\Delta\)¹⁴C-CO₂ = -396.7.

Incubation vs. bulk soil \(\Delta\)¹⁴C

Fig. 10. \(\Delta\)¹⁴C of 2019 bulk soil incubations and corresponding bulk soil

Caption: \(\Delta\)¹⁴C of bulk soil and respired CO₂ by depth for each site in 2019. Error bars show one standard deviation for bulk soil, points show mean of three replicate profiles for bulk soils and single observations for respired CO₂ (still waiting for data from laboratory duplicates).

Fig. 10. Regression of 2019 bulk soil incubations and corresponding bulk soil \(\Delta\)¹⁴C

Caption: Regressions of \(\Delta\)¹⁴C of bulk soil and respired CO₂ by depth for each site in 2019. Error bars show one standard deviation for bulk soil, points show mean of three replicate profiles for bulk soils and single observations for respired CO₂ (still waiting for data from laboratory duplicates).

Time series: \(\Delta\)¹⁴C by depth (as measured)

Fig. 10. Time series of \(\Delta\)¹⁴C by depth, as measured

Caption: Points show mean of three profile replicates for 2001 and 2019 samples. Error bars show ± 1 standard deviation of the mean (only a single profile was analyzed in 2009).

Spline fitting

Fit splines to 2009 data to predict at 2001 depths

Soils collected in both the 2001 and 2009 sampling campaigns were sampled by horizon, but the depth intervals differed between the two sampling years. In 2009, full profiles were excavated for each site, as opposed to the shorter profiles collected in 2001 from the GR and AN sites. Radiocarbon was measured on all three replicate profiles at each site for the 2001 samples, but only for one of the replicate profiles at each site in 2009, e.g. ANpp rep2, etc.

In order to compare the radiocarbon profiles between 2001, 2009, and 2019 we interpolated both radiocarbon and carbon stock data at 1 cm intervals for each site in the 2009 and 2019 datasets, and then summed the carbon-stock-weighted radiocarbon values for each depth interval sampled in 2001. A monotonic cubic spline fit was used for the carbon stock interpolation (Wendt and Hauser 2013), and a mass-preserving spline was used to fit the radiocarbon data (Bishop, T.F.A., McBratney, A.B., Laslett, G.M., (1999) Modelling soil attribute depth functions with equal-area quadratic smoothing splines. Geoderma, 91(1-2): 27-45).

Fig. 10. Time series of \(\Delta\)¹⁴C by depth (2001, 2009, 2019)

Caption: Points for 2001 samples show the mean \(\Delta\)¹⁴C values at the measured depths. Points for 2009 and 2019 samples are spline-fitted estimates of \(\Delta\)¹⁴C predicted for the same depth intervals as measured in 2001. Error bars show ± 1 standard deviation of the mean of three replicate profiles for 2001 and 2019 samples (only a single profile was analyzed in 2009).

Fig. 10. Time series of \(\Delta\)¹⁴C by depth (splined to 2019 depths)

Caption: Points for 2019 samples show the mean \(\Delta\)¹⁴C values at the measured depths. Points for 2001 and 2009 samples are spline-fitted estimates of \(\Delta\)¹⁴C predicted for the same depth intervals as measured in 2019. Error bars show ± 1 standard deviation of the mean of three replicate profiles for 2001 and 2019 samples (only a single profile was analyzed in 2009). NB: Only two depth intervals were measured at the cool and cold andesite sites (max depth of 27 and 28 cm, respectively), so linear extrapolation (using the slope of the last 1cm spline-fitted depth increment) was used to extend the profiles to 30 cm.

Initial modeling

The goal of this modeling exercise is to see how parent material and climate/ecosystem affect estimates of soil carbon ages and transit times. Bulk soil ¹⁴C observations from 2001, 2009, and 2019 will be used to constrain the carbon models, as well as observations of ¹⁴C-CO₂ from laboratory soil incubations of soils collected in 2001 and 2019. Previous work has indicated that the carbon stocks at these sites is likely at equilibrium, so we will apply the steady-state assumption to the modeling.

Two-pool models

One pool models have been shown repeatedly to be inadequate for describing soil carbon dynamics. However, as simple models are easier to constrain, we will start with a two-pool parallel and two-series models, as these are the simplest model system beyond the single pool approach.

The two-pool parallel model requires the following parameters: * decomposition constants for each pool (k₁, k₂) * input partitioning coefficient (\(\gamma\)) * steady-state carbon stocks (C) * inputs (I) * initial values of ¹⁴C

The two-pool series model requires the following parameters: * decomposition constants for each pool (k₁, k₂) * transfer coefficient (\(\alpha\)) * steady-state carbon stocks (C) * inputs (I) * initial values of ¹⁴C

Decomposition rates (k) are related to the amount of ¹⁴C in a pre-bomb system (fraction modern, F) at steady-state by the following equations (cf. Schuur, Druffle, and Trumbore, 2016): >Eq. 1

\[F = \frac{k}{k + \lambda}\] >Eq. 2

\[k = \frac{\lambda \cdot F}{1 - F}\] >where \(\lambda\) is the radioactive decay constant (1/8267).

As the decomposition rates will vary, the initial ¹⁴C content can be determined dynamically with equation 1.

Carbon stocks are known, while inputs will be estimated and are related to the steady-state conditions by the following equation: >Eq. 3

\[I = (k_{1} \cdot C_{1}) + (k_{2} \cdot C_{2})\] >where C₁ and C₂ are the carbon stocks of the two model pools.

Both stocks and inputs can be scaled to the known value of the total carbon pool once the steady-state parameters (k₁, k₂, and \(\gamma\) or \(\alpha\)) have been determined. Pool sizes are a function of the inputs and input partitioning coefficient at steady-state.

A Monte-Carlo Markov chain approach will be used for parameter estimation in combination with an initial optimization algorithm to determine the best set of initial parameters.

Workflow

Initial model fitting was performed for both model structures using generous parameter ranges [0, 1] for all three parameters (k₁, k₂, \(\gamma\) or \(\alpha\)). The initial parameter set was found by fitting the models by eye, followed by optimization with the function “modFit” (R package FME), using the Nelder-Mead algorithm. The best set of parameters found by modFit was then used as the input to a Monte Carlo Markov Chain (MCMC), using the function “modMCMC” (R package FME). The number of iterations for the MCMC optimization was set at 5000 intially, with delayed rejection employed to increase efficiency.

The sum of the mean squared error for the best parameter set was slightly lower for the parallel structure than for the series structure. Additionally, the overall mean error of the residuals was also lower for the parallel structure, moderately so for the bulk C observations but substantially so for the respiration observations (in andesite and granite soils in particular).

However, these initial fits yielded unrealistic parameter estimates for multiple sites, particularly at the lower depths. Additionally, the modFit output showed very high correlation between the parameters for both model structures (slightly higher for the two-pool series model).

[[1]]

[[2]]

[[3]]

[[4]]

[[5]]

[[6]]

[[7]]

[[8]]

[[9]]

[[1]]

[[2]]

[[3]]

[[4]]

[[5]]

[[6]]

[[7]]

[[8]]

[[9]]

Bayesian parameter estimation (MCMC)

This section needs some work.